Here is a rule set that will do the job:
[ a*(b + c) := a*b + a*c,
opt(a) O(x^n) + opt(b) O(x^m) := O(x^n) :: n <= m
:: constant(a) :: constant(b),
opt(a) O(x^n) + opt(b) x^m := O(x^n) :: n <= m
:: constant(a) :: constant(b),
a O(x^n) := O(x^n) :: constant(a),
x^opt(m) O(x^n) := O(x^(n+m)),
O(x^n) O(x^m) := O(x^(n+m)) ]
If we really want the + and * keys to
operate naturally on power series, we should put these rules in
EvalRules. For testing purposes, it is better to put
them in a different variable, say, O, first.
The first rule just expands products of sums so that the rest
of the rules can assume they have an expanded-out polynomial to
work with. Note that this rule does not mention
‘O’ at all, so
it will apply to any product-of-sum it encounters—this rule
may surprise you if you put it into EvalRules!
In the second rule, the sum of two O's is changed to the smaller O. The optional constant coefficients are there mostly so that ‘O(x^2) - O(x^3)’ and ‘O(x^3) - O(x^2)’ are handled as well as ‘O(x^2) + O(x^3)’.
The third rule absorbs higher powers of ‘x’ into O's.
The fourth rule says that a constant times a negligible quantity is still negligible. (This rule will also match ‘O(x^3) / 4’, with ‘a = 1/4’.)
The fifth rule rewrites, for example, ‘x^2 O(x^3)’ to ‘O(x^5)’. (It is easy to see that if one of these forms is negligible, the other is, too.) Notice the ‘x^opt(m)’ to pick up terms like ‘x O(x^3)’. Optional powers will match ‘x’ as ‘x^1’ but not 1 as ‘x^0’. This turns out to be exactly what we want here.
The sixth rule is the corresponding rule for products of two O's.
Another way to solve this problem would be to create a new
“data type” that represents truncated power series.
We might represent these as function calls
‘series(coefs,
x)’ where
coefs is a vector of coefficients for
‘x^0’,
‘x^1’,
‘x^2’, and so
on. Rules would exist for sums and products of such
series objects, and as an optional convenience could
also know how to combine a series object with a
normal polynomial. (With this, and with a rule that rewrites
‘O(x^n)’ to
the equivalent series form, you could still enter
power series in exactly the same notation as before.) Operations
on such objects would probably be more efficient, although the
objects would be a bit harder to read.
Some other symbolic math programs provide a power series data
type similar to this. Mathematica, for example, has an object
that looks like ‘PowerSeries[x,
x0,
coefs,
nmin,
nmax,
den]’, where
x0 is the point about which the power series is taken
(we've been assuming this was always zero), and nmin,
nmax, and den allow pseudo-power-series
with fractional or negative powers. Also, the
PowerSeries objects have a special display format
that makes them look like ‘2 x^2 +
O(x^4)’ when they are printed out. (See
Compositions, for a
way to do this in Calc, although for something as involved as
this it would probably be better to write the formatting routine
in Lisp.)